r\  o 


'  REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-01-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection 
of  information,  including  suggestions  for  reducing  the  burden  to  Department  of  Defense,  Washington  Headquarters  Services  Directorate  for  Information  Operations  and  Reports 
(0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be 
subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

August  2007  REPRINT 

3.  DATES  COVERED  (From  -  To) 

>- 

a. 

O 


4.  TITLE  AND  SUBTITLE 

Lattice  model  of  fluid  turbulence 


t- 

Q 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 

61102F 


6.  AUTHORS 
Jeffrey  Yepez 
George  Vahala* 
Linda  Vahala** 
Min  Soe# 

Sean  Ziegeler## 


5d.  PROJECT  NUMBER 

2304 


5e.  TASK  NUMBER 

OT 


5f.  WORK  UNIT  NUMBER 

A1 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  /RVBYA 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL- V  S-H  A-TR-2008- 1 059 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 
AFRL/RVBYA 


-M  SPONSOR/MONITOR’S  REPORT 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  distribution  unlimited. 


20080819  231 


13.  SUPPLEMENTARY  NOTES  Reprinted  from:  NAVO  MsRC  Navigator,  Spring  2007,  Naval  Oceanographic  Office  Major  Shared  Resource  LJr., 
Stennis  Space  Center,  MS  39522*Coll.  of  William  and  Mary,  Williamsburg,  VA  23185  **01d  Dominion  Univ.,  Norfolk,  VA  23529  #Rogers 

State  Univ.,  Claremore,  OK  74017  ##High  Performance  Computing  Modernization  Program,  Mississippi  State  University,  MS  38677 


14.  ABSTRACT 

The  kinetic  lattice  gas  is  an  optimal  model  of  a  vast  many-particle  system  (such  as  a  fluid)  with  complicated  particle-particle  interactions  and 
irregular  boundary  conditions.  With  it  the  fluid  dynamicist  can  achieve  higher  nonlinearity  (measured  by  Reynolds  number),  unconditional 
stability  and  accuracy,  with  less  memory  and  processor  time  than  with  other  models  of  turbulence  for  situations  with  tortuous  boundaries.  For  an 
engineer,  it  is  simple  to  code,  runs  perfectly  on  parallel  supercomputers,  and  is  suited  to  a  plethora  of  computational  physics  applications.  As 
demonstrated  here,  it  is  a  competitive  alternative  to  large  eddy  simulations  with  Smagorinsky  sub-grid  closure.  To  theorists  and  experimentalists 
in  quantum  information  science,  its  kinetic  transport  equation  is  a  special  case  of  the  quantum  dynamics,  particularly  governing  a  parallel  array  of 
quantum  processors,  a  type-1 1  quantum  computer  architecture.  Presented  are  turbulent  fluid  simulations  using  the  kinetic  lattice  gas  model  carried 
out  on  the  supercomputer  BABBAGE.  As  an  illustration  of  the  efficiency  of  the  lattice  model,  presented  is  a  discovery  of  a  universal  range  in  the 
morphological  evolution  of  the  laminar-to-turbulent  flow  transition:  the  breaking  subrange. 


15.  SUBJECT  TERMS 

Kinetic  lattice  gas  model 
Type  II  quantum  computing 

Entropic  lattice  Boltzmann  equation 

Fluid  turbulence 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

C.  THIS  PAGE 

ABSTRACT 

OF 

PAGES 

Jeffrey  Yepez 

UNCL 

UNCL 

UNCL 

UNL 

19B.  TELEPHONE  NUMBER  (Include  area  code) 

Standard  Form  298  (Rev  8/98) 

Prescribed  by  ANSI  Std.  Z39  18 


DTIC  COPY 


AFRL-RV-HA-TR-2008-1059 

Navigator  2007  -  Centerfold  -  NAVO  MSRC  CAP  II  -  Version  1.2.3.3  -  Cleared  for  public  release:  VS07.005-129  Case  No.  ESC  07-0689. 

Lattice  model  of  fluid  turbulence 

Jeffrey  Yepez1,  George  Vahala2,  Linda  Vahala3,  Min  Soe4,  Sean  Ziegeler5 
1  Air  Force  Research  Laboratory,  Hanscom  Air  Force  Base,  MA  01731 
2  Department  of  Physics ,  William  &  Mary ,  Williamsburg,  VA  23185 
3  Department  of  Electrical  &  Computer  Engineering,  Old  Dominion  University,  Norfolk,  VA  23529 
4  Department  of  Mathematics  and  Science ,  Rogers  State  University,  Claremore,  OK  74017 
5  High  Performance  Computing  Modernization  Program,  Mississippi  State  University,  MS  38677 

The  kinetic  lattice  gas  is  an  optimal  model  of  a  vast  many-particle  system  (such  as  a  fluid)  with 
complicated  particle-particle  interactions  and  irregular  boundary  conditions.  With  it  the  fluid  dy- 
namicist  can  achieve  higher  nonlinearity  (measured  by  Reynolds  number),  unconditional  stability 
and  accuracy,  with  less  memory  and  processor  time  than  with  other  models  of  turbulence  for  sit¬ 
uations  with  tortuous  boundaries.  For  an  engineer,  it  is  simple  to  code,  runs  perfectly  on  parallel 
supercomputers,  and  is  suited  to  a  plethora  of  computational  physics  applications.  As  demonstrated 
here,  it  is  a  competitive  alternative  to  large  eddy  simulations  with  Smagorinsky  sub-grid  closure. 

To  theorists  and  experimentalist  in  quantum  information  science,  its  kinetic  transport  equation  is 
a  special  case  of  the  quantum  dynamics,  particularly  governing  a  parallel  array  of  quantum  pro¬ 
cessors,  a  type-II  quantum  computer  architecture.  Presented  are  turbulent  fluid  simulations  using 
the  kinetic  lattice  gas  model  carried  out  on  the  supercomputer  Babbage.  As  an  illustration  of  the 
efficiency  of  the  lattice  model,  presented  is  a  discovery  of  a  universal  range  in  the  morphological 
evolution  of  the  laminar-to-turbulent  flow  transition:  the  breaking  subrange. 

Keywords:  kinetic  lattice  gas  model,  type  II  quantum  computing,  entropic  lattice  Boltzmann  equation,  fluid 
turbulence 


I.  INTRODUCTION 

In  the  defense  community  there  are  a  plethora  of 
neutral-fluid  flows  problems,  particularly  flows  by  regions 
with  non-trivial  spatial  boundaries,  such  around  an  air¬ 
craft  fuselage  or  a  ship  hull,  through  a  jet  engine  com¬ 
partment,  or  over  the  Earth’s  surface.  While  analyti¬ 
cal  solutions  to  such  problems  remain  elusive  and  gen¬ 
erally  intractable,  today  they  are  within  the  reach  of 
the  computational  physicist.  To  tackle  such  problems, 
practitioners  usually  resort  to  direct  numerical  simula¬ 
tion  methods,  yet  here  the  amount  of  computer  memory 
and  processing  time  grows  faster  that  the  number  of  de¬ 
sired  computed  field  points.  Even  in  cases  with  simple 
or  periodic  boundaries,  where  more  efficient  numerical 
representations  are  available  (such  as  a  psuedo  spectral 
approach  pioneered  by  G.I.  Taylor  in  the  1930’s  and  S.A. 
Orszag  et  al.  in  the  1970’s  [1]),  the  scaling  of  computer 
resources  with  grid  size  is  daunting.  This  computational 
complexity  has  translated  into  significant  annual  costs  to 
defense  department  high  performance  computing  offices 
purchasing  large  numbers  of  processing  elements  (typi¬ 
cally  thousands  or  ten  of  thousands),  configured  in  mas¬ 
sive  parallel  computing  arrays.  But  this  cost,  increasing 
year  after  year,  has  been  borne  so  engineers  can  solve  mis¬ 
sion  critical  fluid  problems,  perpetually  requiring  higher 
resolution  grids  and  more  accurate  and  faithful  simula¬ 
tions.  For  example,  high  resolution  flow  simulations  are 
vital  to  aeronautical  engineers  designing  the  shape  of  ad¬ 
vanced  fighter  jets  or  unmanned  aerial  vehicles  or  sub¬ 
marines  to  economize  on  fuel  consumption  and  minimize 
maneuvering  instabilities  and  wakes,  to  propulsion  engi¬ 
neers  designing  nozzle  and  flow  control  orifices  to  max¬ 
imize  thrust,  or  to  meteorologists  trying  to  understand 


intermittent  turbulence  induced  in  the  upper  atmosphere 
under  the  jet  stream  to  maximize  laser  propagation  from 
airborne  platforms.  Yet  to  the  theoretical  physicist,  the 
situation  is  even  more  dire:  the  prediction  of  any  aspects 
of  turbulence  (beyond  Kolmogorov’s  1941  universality 
hypothesis),  using  advanced  statistical  methods  and  per¬ 
turbation  methods,  borrowed  from  triumphant  quantum 
field  theory  and  statistical  mechanics,  remains  the  oldest 
and  most  prominent  of  classical  grand  challenge  prob¬ 
lems,  open  now  for  over  150  years. 

This  dire  situation  arises  because,  even  in  the  macro¬ 
scopic  limit,  strong  correlations  and  feedback  mecha¬ 
nisms  between  large  scale  and  small  scale  flow  struc¬ 
tures,  over  many  decades  of  spatial  separation,  dominate 
the  overall  flow  evolution.  The  clearest  high  level  pic¬ 
ture  capturing  the  essential  physics  of  this  problem,  with 
restrictioned  attention  to  divergence  free  and  low  Mach 
number  flows,  are  the  incompressible  Navier-Stokes  equa¬ 
tions.  The  strong  correlation  between  disparate  scales  is 
captured  by  the  extremely  simple  non-local  convective 
derivative  (the  second  order  nonlinearity  in  the  velocity 
field). 


II.  LAMINAR  TO  TURBULENT  FLOW 
TRANSITION 

The  lattice  model  now  affords  a  deep  insight  into  the 
origin  and  essential  inner  workings  of  free  shear  turbu¬ 
lence.  This  is  a  kinetic  lattice  gas  model ,  the  clearest 
low  level  picture  correctly  capturing  the  essential  physics 
and  hydrodynamics  of  the  problem.  As  the  well  known 
Ising  lattice  gas  model  is  fundamental  to  a  statistical  me¬ 
chanics  understanding  of  the  essential  physics  of  ferro- 
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FIG.  2:  Surfaces  of  constant  enstrophy  |  /dl/|u;|2  where 

u;  =  V  x  u)  illustrating  an  incompressible  fluid's  morphological 
evolution  from  t  =  0  up  to  t  =  7,  000 A f  iterations,  in  time  steps 
of  l,000Af  on  a  cubical  cartesian  grid  of  size  L  =  5 12 Ax.  Surface 
coloring  uses  u  •  u  (red  equals  -1  and  blue  1).  This  3+1  dimen¬ 
sional  turbulent  neutral  fluid  simulation  run  was  on  the  newest 
defense  department  supercomputer,  Babbage,  using  the  entropic 
lattice  Boltzmann  equation  with  15-body  particle-particle  collisions 
(ELB-Q15  model)  computed  at  every  lattice  site  at  each  time  step. 
At  each  site,  local  relaxation  of  the  single-particle  probability  dis¬ 
tribution  a  desired  equilibrium  function,  represented  as  a  low  Mach 
number  polynomial  expansion.  The  initial  flow  is  a  Kida  and 
Murakami  profile  [2]  with  a  super  cell  size  set  to  L0  =  5 12 Ax, 
the  total  grid  size.  So  the  flow  configurations  within  all  8  oc¬ 
tants  of  the  large  grid  are  initially  identical.  The  characteristic 
flow  speed  is  u0  =  0.07071  The  collisional  inversion  parame¬ 
ter  is  set  to  (3  =  0.99592,  corresponding  to  a  kinematic  viscosity 
of  i/0  =  6.8  x  10~4^^-,  for  q  =  2.  The  Reynolds  number  is 
Re  =  =  53, 024.  The  resulting  turbulence  is  not  fully  re¬ 

solved  down  to  the  dissipation  scale,  which  in  the  model  is  the  cell 
size  Ax.  Do  do  this,  set  L  =  Re^  ~  2, 078 Ax.  So  the  flow  is  under 
resolved  by  about  a  factor  of  4. 


FIG.  1:  Supercomputer  simulation  of  Navier-Stokes  turbulence. 
Surface  of  constant  vorticity  in  turbulent  flow,  showing  entangled 
vortex  tubes  (top)  and  anisotropic  flow  in  the  breaking  subrange 
(bottom).  The  dot  product  of  the  velocity  and  vorticity  fields  are 
displayed  in  the  red-blue  color  coding.  The  bottom  image  is  a 
zoomed  view  of  at  t  =  5KAt  time  step  shown  if  Fig.  2,  where  the 
vertical  white  line  is  one  edge  of  the  cubical  simulation  grid.  Long 
vortex  tendrils  are  easily  seen.  Fig.  1  shows  the  fastest  executing 
code  for  computational  fluid  dynamics  simulations  of  turbulence 
to  date.  This  was  computed  on  the  newest  defense  department 
supercomputer,  Babbage,  located  at  the  Naval  Oceanographic  Of¬ 
fice  MSRC  at  the  Stennis  Space  Center  in  Mississippi,  through  a 
Capability  Application  Program  (CAP)  II  grant.  Also  tested  were 
new  sub-grid  models  of  turbulence,  using  lattice  Boltzmann  equa¬ 
tion  techniques,  an  entropic  method  and  a  Smagorinsky  closure 
method.  These  codes  are  close  cousins  to  quantum  algorithms  for 
aerodynamics.  Our  CAP  II  simulations  used  over  one  million  hours, 
using  a  dedicated  block  of  thousands  of  high  performance  proces¬ 
sors  over  the  period  of  four  weeks,  generating  terabytes  of  data  per 
day.  These  large-scale  simulations  offer  a  better  understanding  of 
the  morphological  evolution  and  structural  development  of  turbu¬ 
lence  in  fluids,  but  they  also  give  a  preview  of  the  kind  of  numerical 
output  to  be  available  from  future  quantum  computers. 
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magnetism  and  the  order-disorder  phase  transition,  the 
kinetic  lattice  gas  model  is  fundamental  to  a  dynamical 
mechanics  understanding  of  fluids  and  the  laminar-to- 
turbulent  flow  transition.  See  Appendix  A  for  a  brief 
mathematical  overview  of  the  model. 


FIG.  3:  Plot  of  enstrophy  versus  time  (smooth  black  curve)  show¬ 
ing  three  stages  in  the  morphological  evolution:  (1)  vortex  stretch¬ 
ing  range  ( t  <  3, 200A*b  (2)  breaking  subrange  (3,200  <  t  < 
9,000At),  and  (3)  inertial  subrange  (t  >  9,000At).  The  enstrophy 
is  normalized  so  that  at  t  =  0  it  is  unity.  The  isovalues  used  to 
visualize  the  8  images  in  Fig.  2  are  shown  (black  squares).  Stage 
I:  The  initial  exponential  increase  (blue  curve)  of  enstrophy  des¬ 
ignates  the  initial  vortex  stretching  range  with  characteristic  lam¬ 
inar  flow.  There  is  excellent  agreement  between  the  analytical  fit 
(blue  curve)  and  the  enstrophy  data  (black  curve).  STAGE  2:  The 
time  derivative  of  the  enstrophy  curve  (jagged  black  curve)  is  also 
plotted.  The  time  period  of  generally  negative  slope  of  the  enstro¬ 
phy  derivative  (gray  shaded  region)  is  here  termed  the  breaking 
subrange ,  where  large  scale  anisotropic  ordering  of  turbulence  oc¬ 
curs  and  intermittently  breaks  down  over  time.  The  first  major 
breaking  point  occurs  at  about  t  =  3,  200Af  (vertical  red  line)  and 
subsequent  intermediate  avalanches  occur  at  about  t  =  5,  lOOAt 
and  t  =  6, 750A t  (thin  vertical  red  lines),  respectively.  Stage 
3:  The  final  exponential  decrease  of  enstrophy  (red  curve)  desig¬ 
nates  the  inertial  subrange  where  the  homogeneous  and  isotropic 
'‘small  scale”  turbulent  flow  morphology,  with  characteristic  entan¬ 
gled  vortex  tubes,  is  organized  in  a  spatially  self-similar  way.  Here 
the  energy  spectral  density  obeys  the  Kolmogorov  universality  hy¬ 
pothesis,  the  famous  k~  5  power  law  for  energy  cascade  downward 
to  smaller  scales.  The  onset  of  the  inertial  subrange  occurs  close  to 
t  =  9,  000 A t  (vertical  blue  line).  Here  the  velocity  probability  dis¬ 
tribution  functions,  for  each  component,  are  Gaussian  (top  inset) 
and  the  vorticity  probability  distribution  function  approaches  an 
exponential  (bottom  inset).  There  is  excellent  agreement  between 
the  analytical  fit  (blue  curve)  and  the  enstrophy  data  (black  curve). 
There  exists  a  fourth  stage  of  the  morphology  of  turbulence  at  late 
times  (t  >  14,000A£),  not  shown  here,  called  the  ihscous  subrange. 

A  new  universal  feature  of  the  laminar- to- turbulent 
transition  becomes  clear  in  high  Reynolds  number  simu¬ 
lations.  Following  an  initial  period  of  laminar  flow  with 
vortex  sheet  stretching,  and  preceding  the  final  inertial 
subrange  period  of  isotropic  and  homogeneous  turbulent 
flow  with  self-similar  vortex  tubes,  there  exists  a  well- 
defined  interim  period,  here  termed  the  breaking  sub¬ 
range.  The  hypothesis  is  that  this  subrange  manifests 
self-organized  criticality.  The  breaking  subrange  is  dom¬ 


inated  by  anisotropic  and  noil-homogeneous  turbulent 
flow.  Avalanches  occur  intermittently,  where  large  coher¬ 
ent  spatial  structures  grow,  become  unstable  under  max¬ 
imal  shear,  and  subsequently  break  into  isotropic  and 
homogeneous  turbulence.  These  avalanches  occur  pro¬ 
gressively  in  time  across  the  entire  space.  See  the  bottom 
image  of  Fig.  1  showing  the  kind  of  anisotropic  turbulent 
flow  that  occurs  during  an  avalanche  event,  with  young 
vortex  tendrils.  The  time  rate  of  change  of  the  enstrophy, 
<  has  a  generally  negatively  sloped  avalanche  cascade, 
with  marked  peaks,  indicating  the  successive  breaking 
points,  as  shown  in  Fig.  3. 

By  carefully  comparing  the  renderings  in  Fig.  2  to  the 
enstrophy  plots  in  in  Fig.  3,  it  is  possible  to  see  three  dis¬ 
tinct  morphological  stages  of  the  flow’:  vortex  stretching 
range,  breaking  subrange,  and  the  inertial  subrange.  The 
morphological  evolution  transitions  from  (1)  large  vortex 
sheets  to  (2)  convoluted  vortex  sheets  with  virgin  vortex 
tendrils  to  (3)  small  entangled  vortex  tubes. 


III.  Q  VERSUS  S  MODELS:  CLOSURE  OF 
SUB-GRID  EFFECTS 

The  continuum  hydrodynamical  equations  are  projec¬ 
tions  of  the  entropic  lattice  Boltzmann  (ELB)  equation, 
a  projection  down  from  the  Q  x  LD-dimensional  kinetic- 
phase  space  on  to  a  4  x  LD-dimensional  hydrodynamic 
null  space  l.  This  projection  recovers  the  Navier-Stokes 
equations  in  the  Chapman-Enskog  limit.  This  has  an 
important  consequence:  there  exist  many  “qubit”  mod¬ 
els  (or  Q  models  for  short)  with  a  different  local  sten¬ 
cils  ( i.e .  lattice  vectors  sets  with  different  finite  point 
group  symmetries  and  coordination  number),  which  will 
also  recover  the  Navier-Stokes  equations  asymptotically 
(continuous  rotational  symmetry). 

ELB  is  ideal  for  large  eddy  simulation  (LES)  closures 
since  in  LES  one  typically  deals  with  mean  strain  rates 
for  modeling  the  eddy  viscosity.  These  nonlocal  fluid 
calculations  are  immediately  recovered  from  simple  local 
moments  in  ELB. 

At  this  stage,  the  ELB  runs  are  a  validation  of  the 
method  and  this  is  important  because  ELB  is  a  crucial 
precursor  to  a  viable  quantum  lattice  Boltzmann  method 
[3].  The  good  numerical  agreement  between  the  LES-LB 
(lattice  Boltzmann  equation  with  sub-grid  Smagorinsky 
closure  here  referred  to  as  an  S  model)  and  the  basic  ELB 
Q  models  is  encouraging.  Now  fully  convinced  of  the  va¬ 
lidity  of  ELB,  work  to  speed  it  up  for  present  day  super¬ 
computer  implementations  is  underway.  But  the  most 


A  subset  of  all  the  eigenvectors  in  the  Q-dimensioual  kinetic 
space,  have  zero  eigenvalues.  These  eigenvectors  span  the  hy¬ 
drodynamic  null  space.  It  is  4  dimensional  because  the  kinetic 
lattice  gas  model  conserves  probability  (mass)  and  probability 
flux  along  the  spatially  orthogonal  directions  (three  components 
of  the  momentum  vector). 
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promising  opportunity  is  to  simply  build  a  type-II  quan¬ 
tum  computer  [4,  5],  which  could  far  outstrip  any  classi¬ 
cal  supercomputers,  for  turbulence  fluid  simulations.  In 
the  fullness  of  time,  the  ELB  will  outpace  the  LES-LB, 
even  without  quantum  computers. 

There  are  a  few  reasons  for  this  view.  First,  the  kinetic 
lattice  gases  obey  detail  balance  while  sub-grid  closure 
methods,  such  as  the  LES-LB,  do  not.  And  another  ad¬ 
vantage  of  the  ELB  over  LES-LB  is  that  ELB  obeys  the 
second  law  of  thermodynamics  while  LES-LB  and  other 
LES  methods  do  not  necessarily  obey  the  second  law. 
Third,  in  the  LES-LB  the  strain  tensor  must  be  com¬ 
puted  at  every  site.  With  no-slip  boundaries,  computing 
the  strain  tensor  becomes  problematic.  In  contrast,  ELB 
is  purely  local,  so  grid  sites  near  boundaries  are  handled 
as  easily  as  sites  far  away  from  boundaries.  All  these  are 
important  differences  when  the  model  is  used  for  practi¬ 
cal  engineering  grade  applications. 


IV.  TIMINGS  AND  SCALING 

Turbulent  dynamics  are  easily  solved  since  the  under¬ 
lying  kinetic  equation  (Al)  has  only  local  algebraic  non- 
linearities  in  the  macroscopic  variables  and  simple  linear 
advection.  At  this  mesoscopic  level  there  are  various  ki¬ 
netic  lattices  (Q=15,  19,  or  27)  with  different  lattice  vec¬ 
tors  on  a  cubic  lattice,  which  model  the  Navier-Stokes 
equation  to  leading  order  in  the  Chapman-Enskog  per¬ 
turbative  asymptotics. 

With  the  CAP  data  obtained  on  Babbage,  the  effects 
of  the  underlying  lattice  symmetry  on  the  fluid  turbu¬ 
lence  statistics  (through  autocorrelation  tensors  of  veloc¬ 
ity,  vorticity,  pdfs  of  vorticity,  and  the  like)  can  be  deter¬ 
mined,  but  there  is  not  have  sufficient  space  to  present 
details  here.  The  Q15  model  seems  to  be  the  most  effi¬ 
cient  model.  An  example  output  of  this  model  is  shown 
in  Fig.  2.  Even  on  a  relatively  modest  size  5123  grid,  we 
can  achieve  such  a  high  degree  of  resolved  nonlinearity 
(Re=26,512)  that  the  consequent  isotropic  turbulence  at 
the  onset  of  the  inertial  subrange  (at  about  t  ~  9,  OOOAf ) 
outstrips  the  ability  to  visualize  the  myriad  vortex  tubes. 
Fig.  2  is  just  a  test  case.  The  sustained  floating  point 
per  seconds  (MFLOPS/PE)  we  achieved  are  the  best  of 
all  scientific  codes  run,  for  example  on  the  Earth  Simu¬ 
lator,  attaining  over  67%  of  peak  performance  on  4,800 
processors.  Achieving  the  world’s  highest  Reynolds  num¬ 
ber  in  the  field  of  computational  fluid  dynamics  should 
also  occur  in  the  near  future.  The  current  tested  code 
on  Babbage,  on  1.6003  grid  with  2,048  processing  ele¬ 
ments,  can  achieve  a  Reynold’s  number  of  Re=565,667. 
Modeling  atmospheric  scale  turbulence,  in  the  range  of 
Re  ~  106  is  possible  today,  for  the  first  time  in  the  half 
century  long  history  of  numerical  digital  computers  ap¬ 
plied  to  aerodynamics. 

Some  runs  with  these  models-in  particular  the  Q15, 
Q19,  Q27  and  S27  models-have  been  completed.  The 
advantages  of  these  lower  Q  models  are  reduced  wall- 
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FIG.  4:  TFlops/sec  scaling  of  ELB-Q27  code  on  Babbage  with 
number  of  CPUs. 


clock  times  with  less  memory  demands.  The  nonlinear 
convective  derivatives  of  in  the  Navier-Stokes  equation 
are  recovered  from  purely  local  moments  of  kinetic  space 
distribution  function.  This  is  the  basic  reason  why  ELB 
scales  so  well  with  PEs:  the  algorithm  consists  of  simple 
local  computations  and  streaming  of  information  only  to 
near  neighbor  grid  sites. 


No.  PEs 

GRID 

MODEL 

WALLCLOCK  (s) 

GFlops/s  per  PE 

2912 

19503 

ELB-Q27 

7,554.7 

2.17 

2912 

19503 

ELB-Q19 

5.602.7 

2.24 

2912 

19503 

ELB-Q15 

4,798.4 

2.05 

2912 

19503 

LES-LB-S27 

4,451.2 

1.05 

TABLE  I:  The  gigaflops  per  second  per  processor  element  for 
2912  CPU  runs  on  a  1952  x  1946  x  1950  grid  for  four  lattice 
Boltzmann  codes  variants.  The  wallclock  time  is  for  2,000Af 
(lattice  time  steps).  A  full  turbulence  simulation  takes  about 
54,000  time  steps. 

During  Phase  I,  investigating  the  scaling  properties 
of  the  Q27  code,  over  6.3  tera  flops  per  second  on  the 
full  2912  processor  elements  available  on  Babbage  was 
achieved,  see  Fig.  4. 

The  LES-LB-S27  code,  which  no  longer  needs  the  so¬ 
lution  of  an  entropy  constraint  equation,  has  also  been 
tested  in  Phase  I.  It  is  less  computationally  intensive  (due 
to  the  avoidance  of  log-calls  and  the  need  for  a  Newton- 
Raphson  root  finder  at  each  spatial  node  and  time  itera¬ 
tion)  and  shorter  wallclock  time  than  the  ELB-Q27  code, 
see  Table  I. 


V.  QUANTUM  INFORMATION  PROCESSING 

Quantum  information  processing  (QIP)  and  quantum 
communications  will  be  integral  to  the  21st  century.  For 
many  years,  QIP  has  been  included  in  the  Developing 


FIG.  5:  Shown  here  is  a  newsworthy  achievement  of  a  first  quan¬ 
tum  information  processor  circuit  to  embody  a  qubit  (left),  a  basic 
device  for  storing  quantum  information,  designed  by  Terry  Orlando 
of  MIT  and  build  at  MIT  Lincoln  Laboratory  at  Hanscom  AFB  in 
2000  [6,  7].  The  circuit  must  be  placed  in  a  dilution  refrigera¬ 
tor.  Air  Force  Office  of  Scientific  Research  (AFOSR)  supported 
novel  quantum  computing  technology  based  on  superconductive 
electronics  under  the  Quantum  Computation  for  Physical  Model- 
ing  (QCPM)  theme,  and  this  research  has  recently  found  follow-on 
use.  A  superconducting  wire  loop  with  multiple  Josephson  junc¬ 
tions  forms  a  qubit  and  such  qubits  are  coupled  together  to  make 
quantum  logic  gates.  AFOSR  funded  this  new  solid-state  tech¬ 
nology,  fabricated  at  the  superconductive  electronics  foundry  at 
MIT  Lincoln  laboratory.  This  helped  establish  the  basic  fabrica¬ 
tion  techniques  to  build  scalable  quantum  computers  and  mapped 
the  quantum  control  methodology,  proven  with  NMR  spectroscopy 
[4,  5],  onto  the  field  of  superconductive  electronics  for  quantum  in¬ 
formation  processing,  mapping  pulse  protocols  to  allow  qubit-qubit 
logical  operations,  constituting  a  basic  2-qubit  quantum  processor. 
Going  from  a  2-qubit  processor  to  a  16-qubit  processor  necessarily 
entails  a  significant  applied  research  effort  recently  announced  by 
D-Wave,  a  Canadian  start-up  company.  D- Wave’s  16-qubit  proces¬ 
sor,  fabricated  by  NASA,  is  shown  (right). 


Science  and  Technologies  List,  in  the  section  of  critical 
information  systems  technology.  There  now  exists  the 
possibility  of  the  development  of  large  quantum  computer 
arrays,  potentially  outstripping  any  supercomputers  now 
used  for  defense  department  computations. 

On  Babbage  Q27,  Q19,  Q15  lattice  Boltzmann  codes 
were  tested,  like  the  one  shown  in  Figure  1.  Using  2,048 
processors  it  took  two  days  to  complete  a  single  job.  A 
1.0243  grid  takes  about  44  hours  for  Q27  on  512  proces¬ 
sors,  while  the  corresponding  run  for  Q15  takes  about 
28  hours.  The  cost  of  the  largest  supercomputer  parallel 
arrays  annually  adds  to  a  significant  fraction  of  a  billion 
dollars  for  new  government-owned  systems  in  the  United 
States  ( e.g .  19,000  processor  Franklin  at  DOE/NERSC 
cost  about  $50  million  and  occupies  the  space  of  a  gym¬ 
nasium).  Remarkably,  exploiting  quantum  mechanical 
complexity  new  quantum  device  technology  can  be  used 
to  efficiently  compute  the  collision  operator  of  the  Q15 
lattice  Boltzmann  code,  the  basic  engine  of  the  code. 

Therefore,  a  relatively  low  cost  and  small  type-II  quan¬ 
tum  computer  with  a  few'  thousand  16-qubit  quantum 
processor  chips,  perhaps  cooled  with  dilution  refrigera¬ 
tion,  could  handle  the  same  supercomputing  job.  It  could 
cost  a  few  million  dollars  to  assemble,  a  couple  orders  of 
magnitude  less  expensive  than  classical  digital  electronics 


based  supercomputers,  such  as  Babbage,  and  physically 
smaller  by  many  orders  of  magnitude  as  well.  Further¬ 
more,  future  quantum  processor  arrays  with  more  qubits 
per  node,  if  available  in  the  near  future,  could  outstrip 
any  traditional  parallel  supercomputer  purchased  under 
the  department  of  defense  high  performance  computing 
modernization  program. 

The  kinetic  lattice  gas  model  has  proven  to  be  a  state- 
of-the-art  tool  for  understanding  the  morphological  evo¬ 
lution  of  turbulence. 
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Appendix  A:  Mathematical  formulation 

In  the  lattice  model,  dynamics  is  projected  into  a  dis¬ 
crete  kinetic  phase  space.  The  logical  “1”  state  (the  ex¬ 
cited  state)  of  a  qubit  | q)  associated  with  the  spacetime 
point  (x,  t)  encodes  the  probability  fq  of  the  existence  of 
a  particle  at  that  point  moving  with  velocity  cq  = 
where  Axq  are  lattice  vectors,  for  q  =  1, 2, . . . ,  Q.  A 
fundamental  property  of  the  lattice  model  is  that  parti¬ 
cle  motions  in  momentum  space  and  position  space  oc¬ 
cur  independently  [8].  Particle  momentum  and  position 
space  motions  are  generated  by  the  combination  of  an 
engineered  qubit-qubit  interaction  Hamiltonian  H'  and 
a  free  Hamiltonian  —ih'52qcq  •  V,  respectively. 

All  the  particle-particle  interactions,  the  2-bodv  up  to 
and  including  all  the  (Q-2)-body  interactions,  generated 
by  H'  are  mapped  to  a  local  collision  function  H '  ttq 
that  depends  on  all  the  fq  s  at  the  lattice  site  [3]. 

In  the  type-II  quantum  computing  case,  quantum  en¬ 
tanglement  is  localized  among  qubits  associated  with  the 
same  (x,£)  [9],  so: 

fq(x,t)  =  fq(x,t)  +  nq(fi,f2,...fQ)  (Ala) 
fq(x,t)  =  fq(x-  -  At).  (Alb) 

where  fq  and  /,'  are  called  the  incoming  and  outgoing 
probabilities,  respectively.  In  the  classical  limit,  there 


V 


exists  a  fundamental  entropy  function 

Q 

Wi,---  =  (A2) 

9=1 

where  the  are  self-consistently  determined  positive 
weights  (Y2q  lq  =  1)-  in  (Ala)  is  determined  by  the 
constant  entropy  condition: 

^(/! . /q)  =  ^(/i,..,M  (A3) 

It  is  (Al)  through  (A3)  that  constitutes  the  lattice  model 
of  fluid  turbulence  suited  for  parallel  supercomputers. 

In  the  Q-dimensional  kinetic  space,  the  equilibrium 
distribution  /£q,  taken  as  a  vector,  is  the  bisector  of  the 
difference  of  the  incoming  and  outgoing  kinetic  vectors: 

4-C  +  (i-^)  «..<**>> 

when  limr/_+o  ol(3  =  2  and  where  /£q  is  analytically  deter¬ 
mined  by  extremizing  (A3)  subject  to  the  two  constraints 
of  conversation  of  probability  and  probability  flux.  Elim¬ 
inating  Clq  from  (A4)  yields  an  operative  collision  equa¬ 
tion  simpler  than  (Ala) 

/'=/,  + a/3  (/*«-/,)•  (A5) 

Finally,  eliminating  /r'  in  (Alb)  and  (A5),  yields  the  lat¬ 
tice  Boltzmann  equation 

fq(x  +  Axq,t  +  At)  =  fq(x.t)+a(3  [f?(x,t)  -  fq{x,t)]  , 

(A6) 


in  the  BGK  collisional  limit  [10].  Since  a/3  and  /£q  are 
determined  by  (A2),  (A6)  is  called  the  entropic  lattice 
Boltzmann  equation  [11]. 

The  kinetic  lattice  gas  model  becomes  equivalent  to 
the  Navier  Stokes  equations: 


dtu  +  u- S/a  =  -VP  +  r/V1 2 3 4 5 6 7a  and  Vu  =  0,  (A7) 


where  u  =  u(x,t)  is  the  vectorial  flow  field,  where  the 
pressure  is  proportional  to  the  field  density  (P  =  pc2a , 
with  spatial  dimension  D  and  sound  speed  c.s  =  ^37), 
and  where  the  shear  viscosity  is  the  measure  of  dissipa¬ 
tion  (a  renormalized  transport  coefficient  for  momentum 
diffusion).  I11  the  limit  when  A r  — ►  0  and  At  — ►  0  and 
the  hydrodynamic  variables  are  cast  as  moments  of  the 
probability  distribution,  (pc,u)  =  £ „{c,cq)fq  [12].  The 
hydrodynamic  variables  are  independently  evaluated  at 
each  spacetime  point  (x,  t).  The  shear  viscosity  is  ana¬ 
lytically  determined,  and  the  result  is 


1  =  -  5)  ,  (AS) 


so  the  model  approaches  the  inviscid  limit  where  r/  — >  0 
as  a/3  — ►  2  [13]. 
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